R Notebook

# code to load packages
library(tidyverse)
library(tidymodels)
library(knitr)
library(dplyr)
library(ggplot2)
library(patchwork)
library(corrplot)

corrplot 0.94 loaded

#baseR

auto <- read_csv("data/autopod_data.csv")

New names:
• `notes` -> `notes...22`
• `notes` -> `notes...60`
• `` -> `...98`
Rows: 582 Columns: 122
── Column specification ─────────────────────────────────────────────────────────
Delimiter: ","
chr (108): Collector, Higher taxon, Taxon, miss mC, miss mT, miss  man-D, mis...
dbl   (2): radius L, Ulna L
lgl  (12): m2-5 av, hpp2-5 av, hmp2-5 av, hdp2-5 av, notes...60, fpp2-5 av, f...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

#Set hms as numeric 
auto$hm1 <- as.numeric(as.character(auto$hm1))

Warning: NAs introduced by coercion

auto$hm2 <- as.numeric(as.character(auto$hm2))

Warning: NAs introduced by coercion

auto$hm3 <- as.numeric(as.character(auto$hm3))

Warning: NAs introduced by coercion

auto$hm4 <- as.numeric(as.character(auto$hm4))

Warning: NAs introduced by coercion

auto$hm5 <- as.numeric(as.character(auto$hm5))

Warning: NAs introduced by coercion

#Set hpps as numeric 
auto$hpp1 <- as.numeric(as.character(auto$hpp1))

Warning: NAs introduced by coercion

auto$hpp2 <- as.numeric(as.character(auto$hpp2))

Warning: NAs introduced by coercion

auto$hpp3 <- as.numeric(as.character(auto$hpp3))

Warning: NAs introduced by coercion

auto$hpp4 <- as.numeric(as.character(auto$hpp4))

Warning: NAs introduced by coercion

auto$hpp5 <- as.numeric(as.character(auto$hpp5))

Warning: NAs introduced by coercion

#remove NAs
autoclean <- auto %>%
  filter(complete.cases(hm1, hm2, hm3, hm4, hm5, hpp1, hpp2, hpp3, hpp4, hpp5))

#avg hand metatarsals for each higher taxa

hm_avg <- autoclean |>
  group_by(`Higher taxon`) |>
  summarize(average_hm1 = mean(hm1, na.rm = TRUE),
            average_hm2 = mean(hm2, na.rm = TRUE),
            average_hm3 = mean(hm3, na.rm = TRUE),
            average_hm4 = mean(hm4, na.rm = TRUE),
            average_hm5 = mean(hm5, na.rm = TRUE)
            )
print(hm_avg)

#avg proximal phalanges for each higher taxa

hpp_avg <- autoclean |>
  group_by(`Higher taxon`) |>
  summarize(average_hpp1 = mean(hpp1, na.rm = TRUE),
            average_hpp2 = mean(hpp2, na.rm = TRUE),
            average_hpp3 = mean(hpp3, na.rm = TRUE),
            average_hpp4 = mean(hpp4, na.rm = TRUE),
            average_hpp5 = mean(hpp5, na.rm = TRUE)
            )
print(hpp_avg)

#get total digit lengths 

autoclean$total_digit1 <- autoclean$hm1 + autoclean$hpp1
autoclean$total_digit2 <- autoclean$hm2 + autoclean$hpp2
autoclean$total_digit3 <- autoclean$hm3 + autoclean$hpp3
autoclean$total_digit4 <- autoclean$hm4 + autoclean$hpp4
autoclean$total_digit5 <- autoclean$hm5 + autoclean$hpp5

head(autoclean)

NA

#AVG DIGIT LENGTHS
avg_digit_lengths <- hm_avg |>
  full_join(hpp_avg, by = "Higher taxon") |>

  mutate(
      avg_digit1 = hm_avg$average_hm1 + hpp_avg$average_hpp1,
    avg_digit2 = hm_avg$average_hm2 + hpp_avg$average_hpp2,
    avg_digit3 = hm_avg$average_hm3 + hpp_avg$average_hpp3,
    avg_digit4 = hm_avg$average_hm4 + hpp_avg$average_hpp4,
    avg_digit5 = hm_avg$average_hm5 + hpp_avg$average_hpp5
    ) |>
    select(`Higher taxon`, avg_digit1, avg_digit2, avg_digit3, avg_digit4, avg_digit5)

avg_digit_lengths

#EXPLORING COVARIATION
# biplots for comparing digit lengths

#to calc rsquared 
r_squared <- function(x, y) {
  model <- lm(y ~ x)
  summary(model)$r.squared
}

pairs <- list(
  c("total_digit1", "total_digit2"),
  c("total_digit1", "total_digit3"),
  c("total_digit1", "total_digit4"),
  c("total_digit1", "total_digit5"),
  c("total_digit2", "total_digit3"),
  c("total_digit2", "total_digit4"),
  c("total_digit2", "total_digit5"),
  c("total_digit3", "total_digit4"),
  c("total_digit3", "total_digit5"),
  c("total_digit4", "total_digit5")
)

for (pair in pairs) {
  x <- autoclean[[pair[1]]]
  y <- autoclean[[pair[2]]]
  r_squared <- calculate_r_squared(x, y)
  p <- ggplot(autoclean, aes_string(x = pair[1], y = pair[2])) +
    geom_point() +
    labs(x = pair[1], y = pair[2], title = paste("Biplot of", pair[1], "and", pair[2]), subtitle = paste("R^2 =", r_squared)) +
    theme_minimal()
  print(p)
}

#correlation matrix for total digit lengths
correlation_matrix <- cor(autoclean[, c("total_digit1", "total_digit2", "total_digit3", "total_digit4", "total_digit5")])
print(correlation_matrix)

             total_digit1 total_digit2 total_digit3 total_digit4 total_digit5
total_digit1    1.0000000    0.9566183    0.9565192    0.9634344    0.9664044
total_digit2    0.9566183    1.0000000    0.9791062    0.9725032    0.9686648
total_digit3    0.9565192    0.9791062    1.0000000    0.9941848    0.9786560
total_digit4    0.9634344    0.9725032    0.9941848    1.0000000    0.9928006
total_digit5    0.9664044    0.9686648    0.9786560    0.9928006    1.0000000

corrplot(correlation_matrix, method = "color", addCoef.col = "black", tl.col = "black", tl.srt = 45)

# Function to calculate R^2
calc_r_squared <- function(x, y) {
  model <- lm(y ~ x)
  summary(model)$r.squared
}

# Create pairs of digit lengths
pairs <- list(
  c("total_digit1", "total_digit2"),
  c("total_digit1", "total_digit3"),
  c("total_digit1", "total_digit4"),
  c("total_digit1", "total_digit5"),
  c("total_digit2", "total_digit3"),
  c("total_digit2", "total_digit4"),
  c("total_digit2", "total_digit5"),
  c("total_digit3", "total_digit4"),
  c("total_digit3", "total_digit5"),
  c("total_digit4", "total_digit5")
)

# Get unique taxa
unique_taxa <- unique(autoclean$`Higher taxon`)

# for each higher taxa 
for (taxon in unique_taxa) {
  # filter for the current higher taxon
  taxa_data <- autoclean[autoclean$`Higher taxon` == taxon, ]
  
  # going thru each digit pair 
  for (pair in pairs) {
    x <- taxa_data[[pair[1]]]
    y <- taxa_data[[pair[2]]]
    
    # get r^2
    r_squared <- calc_r_squared(x, y)
    
    # actual plot
    p <- ggplot(taxa_data, aes_string(x = pair[1], y = pair[2])) +
      geom_point(color = "#0096FF", alpha = 0.6, size = 3) +
      labs(
        x = pair[1],
        y = pair[2],
        title = paste("Biplot of", pair[1], "and", pair[2], "for", taxon),
        subtitle = paste("R^2 =", round(r_squared, 3))
      ) +
      theme_minimal() +
      theme(
        plot.title = element_text(size = 12, face = "bold", color = "#0096FF"),
        plot.subtitle = element_text(size = 10, color = "black"),
        axis.title.x = element_text(size = 10, face = "bold"),
        axis.title.y = element_text(size = 10, face = "bold"),
        panel.background = element_rect(fill = "#FFF1F3"),  
        panel.grid.major = element_line(color = "white"), 
        panel.grid.minor = element_line(color = "white")
      )
    print(p)
  }
}

unique_taxa <- unique(autoclean$`Higher taxon`)

for (taxon in unique_taxa) {
  taxa_data <- autoclean[autoclean$`Higher taxon` == taxon, ]

  correlation_matrix <- cor(autoclean[, c("total_digit1", "total_digit2", 
                                          "total_digit3", "total_digit4", 
                                          "total_digit5")])
  print(taxon)
  print(correlation_matrix)
  corrplot(correlation_matrix, method = "color", addCoef.col = "black", tl.col = 
             "black", tl.srt = 45,   title=taxon,    mar=c(0,0,1,0) )
}

[1] "Hominoidea"
             total_digit1 total_digit2 total_digit3 total_digit4 total_digit5
total_digit1    1.0000000    0.9566183    0.9565192    0.9634344    0.9664044
total_digit2    0.9566183    1.0000000    0.9791062    0.9725032    0.9686648
total_digit3    0.9565192    0.9791062    1.0000000    0.9941848    0.9786560
total_digit4    0.9634344    0.9725032    0.9941848    1.0000000    0.9928006
total_digit5    0.9664044    0.9686648    0.9786560    0.9928006    1.0000000
[1] "Cercopithecoidea"
             total_digit1 total_digit2 total_digit3 total_digit4 total_digit5
total_digit1    1.0000000    0.9566183    0.9565192    0.9634344    0.9664044
total_digit2    0.9566183    1.0000000    0.9791062    0.9725032    0.9686648
total_digit3    0.9565192    0.9791062    1.0000000    0.9941848    0.9786560
total_digit4    0.9634344    0.9725032    0.9941848    1.0000000    0.9928006
total_digit5    0.9664044    0.9686648    0.9786560    0.9928006    1.0000000

[1] "Dermoptera"
             total_digit1 total_digit2 total_digit3 total_digit4 total_digit5
total_digit1    1.0000000    0.9566183    0.9565192    0.9634344    0.9664044
total_digit2    0.9566183    1.0000000    0.9791062    0.9725032    0.9686648
total_digit3    0.9565192    0.9791062    1.0000000    0.9941848    0.9786560
total_digit4    0.9634344    0.9725032    0.9941848    1.0000000    0.9928006
total_digit5    0.9664044    0.9686648    0.9786560    0.9928006    1.0000000

[1] "Platyrrhini"
             total_digit1 total_digit2 total_digit3 total_digit4 total_digit5
total_digit1    1.0000000    0.9566183    0.9565192    0.9634344    0.9664044
total_digit2    0.9566183    1.0000000    0.9791062    0.9725032    0.9686648
total_digit3    0.9565192    0.9791062    1.0000000    0.9941848    0.9786560
total_digit4    0.9634344    0.9725032    0.9941848    1.0000000    0.9928006
total_digit5    0.9664044    0.9686648    0.9786560    0.9928006    1.0000000

[1] "Scandentia"
             total_digit1 total_digit2 total_digit3 total_digit4 total_digit5
total_digit1    1.0000000    0.9566183    0.9565192    0.9634344    0.9664044
total_digit2    0.9566183    1.0000000    0.9791062    0.9725032    0.9686648
total_digit3    0.9565192    0.9791062    1.0000000    0.9941848    0.9786560
total_digit4    0.9634344    0.9725032    0.9941848    1.0000000    0.9928006
total_digit5    0.9664044    0.9686648    0.9786560    0.9928006    1.0000000

[1] "Strepsirrhini"
             total_digit1 total_digit2 total_digit3 total_digit4 total_digit5
total_digit1    1.0000000    0.9566183    0.9565192    0.9634344    0.9664044
total_digit2    0.9566183    1.0000000    0.9791062    0.9725032    0.9686648
total_digit3    0.9565192    0.9791062    1.0000000    0.9941848    0.9786560
total_digit4    0.9634344    0.9725032    0.9941848    1.0000000    0.9928006
total_digit5    0.9664044    0.9686648    0.9786560    0.9928006    1.0000000

[1] "Tarsiiforms"
             total_digit1 total_digit2 total_digit3 total_digit4 total_digit5
total_digit1    1.0000000    0.9566183    0.9565192    0.9634344    0.9664044
total_digit2    0.9566183    1.0000000    0.9791062    0.9725032    0.9686648
total_digit3    0.9565192    0.9791062    1.0000000    0.9941848    0.9786560
total_digit4    0.9634344    0.9725032    0.9941848    1.0000000    0.9928006
total_digit5    0.9664044    0.9686648    0.9786560    0.9928006    1.0000000